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Abstract 

The solution of the many-body statistical mechanical theory of elasticity formulated by 
James and Guth in the 1940s is presented. The remarkable aspect of the solution is that it gives 
an elastic free energy that is essentially equivalent to that developed by Flory over a period of 
several decades. 

Introduction 

Rubber elasticity is the first bulk property of polymers that yielded to theoretical analysis. The 
identification of the relation between the Gaussian distribution of the end-to-end distance of a 
random walk with the quadratic strain dependence of phenomenological stress-strain theory was 
key to this success. Of course, the underlying physics that makes this connection inevitable and 
viable is that a polymer chain in the bulk amorphous melt phase is unperturbed by intermolecular 
interactions. [HI This fact allows one to realize the force that a chain delivers to the cross-links that 
terminate it, and which tie it to other chains in the three-dimensional space filling random network, 
is determined solely by the chain’s intramolecular potential. While high elasticity theory was 
the first to successfully predict bulk polymeric materials behavior, it remains one of the few, and 
perhaps the only, analytical theory of polymers to do so, computer simulations notwithstanding. 
This is reason enough to justify efforts to improve upon the theory. Given the history of the subject, 
rubber elasticity is one of the first soft materials to admit an atom-based theoretical analysis. 

The extent to which the elastic equation of state is determined by the interaction between net¬ 
work connectivity and chain statistics has been a point of contention from the earliest days of poly¬ 
mer theory. The theory initiated by Kuhn,[|2l elaborated by Wall[|3l|4land Elory and Rehner,[151[6l 
and discussed extensively in treatises, |[T1 171IH HI HOl is constructed by adding together the con¬ 
tributions to the stress from independent chains. This requires the so-called affine assumption - 
the displacement of the ends of an average network chain is congruent to the macroscopic strain. 
This theory is relatively easy to execute, but by treating the chains as independent it incurred the 
criticism of James and Guth. lfTTl fT^ [T3]l In their many-body theory the individual chains obey 
Gaussian statistics, just as in the independent chain theory, but James and Guth emphasized the 
fact that in tying the chains together with cross-links they become an indissoluble whole that must 
be treated as a single entity. This insight carried a heavy price - their many-body formulation 
was too difficult to be convincingly solved. To make progress, James and Guth lfTTI introduced the 
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unphysical notion of fixed junctions (fixed by the external constraints) that are displaced by the 
macroscopic strain. However, the value that one deduces for the modulus depends on the choice 
and number of fixed junctions, lO which effectively ruins any rigor that may be ascribed to the 
many-body theory. 

Improvements to the independent chain theory by Flory and his colleagues over the span of 
about forty years introduced successive improvements to the theory. In treating swelling, Florv lfTSll 
introduced the controversial combinatorial term giving a contribution —(2z//c//) ln(I/) to the en¬ 
tropy. Here v is the number of chains in the network of volume V, f is the functionality of the 
cross-links, and k is Boltzmann’s constant. Using somewhat convoluted reasoning. ifT^ [TTll the 
ratio {rf) /(r^)o was inserted into the theoretical Young’s modulus in the 1950s. The numerator in 
this expression is the value of the mean-square end-to-end distance of the average network chain 
in a reference state and the denominator is the similar quantity for the free unperturbed chain at 
temperature T . By involving chain dimensions in the modulus, Flory and coworkers were able 
to evaluate the temperature dependence of chain dimensions in terms of the stress-temperature co¬ 
efficient. (It will be noted that the unperturbed chain dimensions drop out of the modulus in the 
independent chain theory, so an additional argument is needed to put this term into the modulus.) 
This ratio of dimensions was not featured in Flory’s work ifTdl of the 1970s, where the number 
of chains is replaced by the cycle rank to account for the contribution from cross-links that are 
inserted after the gel point. This appears to have been motivated by the idea that up to the gel 
point the nascent network, which is approximated as an acyclic tree, cannot support an equilibrium 
stress. Given this idea, it is only the cross-links that are inserted after the gel point that contribute 
to the stress. Assuming that all the prepolymer has been incorporated into the tree at the gel point, 
the cross-links that are inserted subsequently can only form cycles in the tree. The cycle rank 
measures the number of chains that are “activated” by these post-gel cross-links. In replacing the 
number of chains by the cycle rank, the modulus decreases to a value close to that advocated by 
James and Guth. 

High elasticity has attracted the attention of theorists far too numerous to mention here. Much 
of this work has been reviewed in articles ifT^ IT9ll and monographs. If9l ITOl This paper is not in¬ 
tended to be a comprehensive review of all the important work that has gone into our current 
understanding of the physics and chemistry of elastomers and gels. It is instead aimed at solving 
the James-Guth many-body theory in the small strain limit. The objective of the paper is to elimi¬ 
nate as many physical assumptions as possible beyond those that comprise the basic model, and to 
make mathematical approximations clear so that appropriate confirmation or improvements might 
be made in the future. The remarkable result of this calculation is that it gives all the terms that 
Flory put in by hand over the years, although the interpretation of the terms is somewhat differ¬ 
ent. There is no essential difference between the James-Guth and Flory-Wall theories! The other 
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aspect of this work to emphasize is that it firmly establishes a baseline against which one may 
make quantitative statements about other models, including those that treat entanglements, e.g., 
tube models. 

Configuration Integral and Potential of Mean Force 

The shape of a soft system is a characteristic of fundamental interest, and unlike a gas or liquid, 
where the volume is fixed by a container, a soft system adopts a macroscopic shape that is deter¬ 
mined by a combination of internal potentials and external constraints that are imposed at the time 
of formation. Once formed, an unconstrained soft system is free to adjust its boundary based only 
on internal potentials. This section formalizes the computation of the probability function that is 
required to specify the geometric information that characterises the size and shape of the body of 
interest. 

The probability that a classical system described by the Hamiltonian H{p, q), and having ther¬ 
mal energy kT, is found in a state {p, g}, where p and q are 3iV dimensional momenta and coordi¬ 
nate vectors, respectively, is proportional to exp[—H{p, q)/kT]. The p and g of a classical system 
are continuously variable, such that the probability P{p, q)dpdq is defined in the 6iV-dimensional 
coordinate patch dpdq by 


P{p, q)dpdq = Z ^ exp[—H{p, q)/kT]dpdq 

where Z normalizes the distribution. Now define a function F(p, g), which in general can be a 
tensor valued function of any degree, with corresponding volume element dF. The probability 
P(F)(iF that the system will be found in a state with value F is given by 

P{F)dF = dF j P{p, q)dpdq/dF, 

with the integration being performed over the space complementary to F, which is the meaning of 

dpdq/dF. 

The change in free energy, AA, accompanying a change in state of the system from Fi to F 2 
is given by 

AA = -kT ln[P(F2)dF2/P(Fi)dFi] = -kT ln[J(F2/Fi)P(F2)/P(Fi)] (1) 

where dF 2 /dFi = J(F 2 /Fi) is the Jacobian determinant of the mapping Fi —)■ F 2 . In most cases 
the Jacobian determinant will not be commensurate with P(F) and will therefore not contribute to 
the thermodynamics. The reversible work done on the system to convert it from state Fi to state 
F 2 is just AA. Because it is only the ratio of probabilities that is important in thermodynamics, the 
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normalization by Z is immaterial and will be dropped where no harm eomes from doing so. 

The elastie body of interest eonsists of a single eovalently bonded gel eomponent, with no 
sol fraetion. For this elassieal problem the momenta integrate trivially and one is left with the 
ealeulation of the unnormalized probability distribution 
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that the system eonsisting of N (implieit) polymer atoms at temperature T has a shape and size 
determined by S = S(g) (to be speeified later). Here V(g) is the potential energy of the system, 
/3 = 1/kT , and dq is the 3iV-dimensional volume element. (At this early stage swelling ean be 
formally aeeommodated with use of a semi-grand ensemble for a polymer-solvent system - the 
system ean be open with respeet to exehange of solvent moleeules at speeified ehemieal poten¬ 
tial. This paper foeuses on the stress-stain relation and not on swelling, so that elaboration is not 
pursued.) 

The “moding-out” operation, dpdq —?• dpdq/dF, is now performed onee more, but on the inside 
of the shape distribution. Let there be p eross-link “atoms”, and eonsider performing the integral 
in eq. (2) in two steps as 



That is, first integrate over the mid-ehain eoordinates between the eross-links while holding the 
latter eoordinates fixed, after whieh the integrations over the eross-link positions are to be exeeuted. 
The first integration generates an aeeeptable approximation to the free energy for the un-eross- 
linked polymer together with a remaining pieee that is determined by the potential of mean foree 
aeting between the eross-links. That is 
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where Ao(S, T) is the free energy of the un-eross-linked polymer. This funetion is presumed to 
depend only on the volume of the elastomer and not on its shape; this will be defined more pre- 
eisely after S is speeified. What remains under the integral is the Boltzmann faetor of the potential 
of mean foree, V{q^), aeting between eross-links. [The Flory-Rehner theory of swelling [|5l O 
assumes this separation of the free energy of the base polymer from that of the elastomer. Sen¬ 
sitive swelling experiments[|^|23|23l suggest that this separation is not strietly valid. A theory 
that eouples eross-link modes of motion with mid-ehain motions would take us far afield of the 
present objeetive, whieh is restrieted to an analysis of the elassieal theory and the presentation of 


4 



techniques for handling soft materials.] 

The potential of mean force, V{q^) , for an elastomer is given in all elementary elasticity 
theories as a product of Boltzmann factors, one for each chain in the network. James and Guth lfTH 
wrote this product as a many-body quadratic potential of the form 

V/kT = Y, nAn - rjf = -fiT{XKX') (4) 

i-j 

where 7 *^ = 3/2(r^)jj is the modulus parameter for the tie-chains, with {r‘^)ij being the mean- 
square displacement, at the implicit temperature T, between junction pairs i — j that are directly 
connected by a single chain. (James and Guth formulated the potential with a distribution of chain 
lengths, as has been done here. The simplified, second version in eq. (4) is written with an average 
7 to emphasize the chain dimension parameter that will carry important information. For a network 
constructed with uniform chains this is the sole required molecular characteristic; for a distribution 
of chain lengths the averaging over the chain length distribution will be addressed later.) The 
Cartesian coordinates of the junction points are written as the 3 x /i matrix X = (xf); I < a < 
3,1 < * < /i, with transpose X'. There may be multiple chains that directly connect two junctions, 
but this is incorporated in the matrix K, which is the Laplacian for the graph [l20l that encodes the 
connectivity of the network. (In previous work the author [l24l named this the Kirchhoff matrix. 
The mathematical literature [l20ll has settled on the name Laplacian. It would also be appropriate 
to call this the Hessian for the network.) Regardless of the name, the construction of the matrix 
has been described in detail in several publications ifTTl l24l |25]| and need not be repeated here. It is 
assumed that the Laplacian has a single zero eigenvalue, which signifies that the system consists 
of a single connected component. 

Brout-Fixman-Edwards (BFE) Averaging 

The next formality turns out to be extremely important for practical evaluation of the stress-strain 
relation: we need to average over frozen disorder. Systems with frozen disorder require a higher 
level average of the free energy than is usually encountered in elementary statistical mechanics. 
This was first described by Brout ll^ in a treatment of order-disorder transitions; Fixman used the 
averaging in the guise of conditional probabilities in work on polypeptides [|271 and polynucleic 
acids : [[28l and finally, Edwards and coworkers lf29ll formalized the averaging in several treatments 
of rubber elasticity. Let {G} be a set of internal constraints. The set {G*} is a particular instantia¬ 
tion of the constraints that characterizes a member of the ensemble of samples of the material, all 
members having been prepared under the action of identical external constraints. The probability 
that a system is observed with this set of constraints is P(T*, {Gj}) oc exp[— 74 (T*, {Gj})/fcT*], 
where A{T^, {G*}) is the free energy of the system that is formed at temperature T* when the con- 
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straints {C} are imposed. While the eonstraints are formed at temperature T*, they remain fixed 
when the system temperature is ehanged to T. The BFE average that we need is 

{A{T)) = y] P(r„ {C,})A(T, {C.}) (5) 

i 

The average {A{T)) depends on additional parameters that are implieit. For the elastieity problem, 
the eure temperature is T*, and the eonstraints are the eross-linkages that are formed at eure. 

The Holonomic Constraint Trick 

The exaet stress-strain relation from eontinuum meehanies lf30l is 

a = 2{pJpi)e[da/d{e'e)]Te' (6) 

where a is the stress tensor, p^ is the mass density at the state of strain speeified by the deformation 
gradient tensor £ (for the unstrained state, e = 1 = unit tensor), e' is the transpose of e, and a is the 
Helmholtz free energy per unit volume in the unstrained state. All of these quantities are defined 
point-wise, and in general will vary from point to point in the medium. While the equation is an 
exaet eontinuum equation, it is untenable from the standpoint of statistieal meehanies. Suppose 
that one has an inhomogeneous medium in whieh the phases are suffieiently finely dispersed that 
the interaetion between neighboring phases is a substantial eontribution to the free energy density. 
How does one evaluate the free energy density for a domain? It is not the diseontinuities per se that 
eause problems - it is the interaetions between neighboring volume elements that defeat evaluation 
of a free energy density. The potential energy of atoms that interaet aeross bounding surfaees of 
volume elements eannot be assigned unequivoeally to the elements on either side of the surfaee. 
This is one of the more glaring examples of the ineompatibility of eontinuum and atomie deserip- 
tions of matter - there are others. Furthermore, for most praetieal applieations the mieroseopie 
strain is neither important nor measurable. The statistieal meehanieal problem is best defined in 
the thermodynamie limit of large systems where the maeroseopie strain is the only eoneern. Onee 
problems at this length seale are solved, the problems of finely divided inhomogeneous media ean 
be approaehed from above, whieh is the usual approaeh in engineering ealeulations with, say, finite 
element methods. 

The many-body theory that James and Guth formulated left them with the diffieult problem of 
relating the moleeular eoordinates to the strain. For a homogeneous erystalline solid, the deforma¬ 
tion of a unit eell eoineides, on average, with the maeroseopie deformation, and this immediately 
relates the eell axes to the maeroseopie deformation. On the other hand, a simple fluid eonforms to 
the shape of its eontainer, so the state of strain is of no eonsequenee. However, eontainers provide 
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non-holonomic constraints on the configuration space of simple fluids: the coordinates are confined 
to a compact domain determined by the container. In usual practice an elastomer is constrained by 
mechanical means over a portion of its surface by a set of clamps, a wheel rim, a road bed, etc., 
which are external constraints that deliver stress. As in the case of fluids, these are non-holonomic 
constraints that impose boundary conditions on the configuration space; for an elastomer these are 
invariably discontinuous boundaries. However, the very complicated constraints that might be en¬ 
countered in a real application should not be solved at the level of statistical mechanics. We have 
to be content to evaluate the equation of state for a simple geometrical shape and leave complex 
geometries to engineering calculations. 

In their 1943 paper lim . James and Guth wrote that “Rubber resembles a gas very strikingly in 
its thermoelastic behavior.” Given this analogy, it may have been natural to think about integrat¬ 
ing over a configuration space with boundary constraints similar to those provided by a container. 
James and Guth introduced holonomic constraints by picking junction points in the network that 
were declared to be fixed by external forces, and which are displaced by the macroscopic strain. 
This artifice enabled them to sidestep difficult integrations. They showed that the average coordi¬ 
nates of the free junctions are linear functions of the coordinates of the fixed junctions; integrations 
over their fluctuating positions have the majority of their support over molecular dimensions. The 
difficult problem of integrating all coordinates over the volume of the elastomer, subject to con¬ 
straints, was thereby eliminated. 

Unfortunately, this construction merely shifted the problem to a different arena. As Flory 
pointed out. lfTdl the value that one obtains for the modulus in this treatment depends on the number 
and location of the fixed points. One cannot make a conclusive theory based on this treatment of 
constraints without additional assumptions. The rigor that is inherent in the many-body theory was 
vitiated by fixing junctions. 

What is needed is a trick to introduce holonomic constraints which convey the dimensional 
information of non-holonomic constraints, but which do not require special treatment for any par¬ 
ticular atomic species. The symmetric gyration tensor 
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provides just what is needed. The continuum form on the right, with integration over the volume 
V of the material, approximates the atomic coordinate form on the left all the better the larger 
(thermodynamic limit) is the system. In both the atomic and continuum forms, the coordinates are 
measured from the center of symmetry of the system. The macroscopic strain is defined by the 


7 


deformation tensor with eomponents via 

^ S = A'SoA (8) 

Here So is a referenee state, the summation eonvention is used, and a matrix representation is 
eonveyed by the boldface equation on the right. It is worthy of note that the equality between 
atomic and continuum evaluations of the gyration tensor in eq. (7) accommodates inhomogeneous 
media as well as homogeneous. The macroscopic strain tensor A defined in eq. (8) is insensitive to 
atomic detail; it is defined by this equation. Note that this is a macroscopic deformation, whereas 
e in eq. (6) is microscopic, i.e., is a function of coordinates (location) in general. 

Eckart Coordinates 

The 3 X yu matrix X = (a;“), 1 < a < 3,1 < z < /r, of cross-link coordinates gives the gyration 
tensor in eq. (7) as S = ^~^XX'. The diagonal components of this tensor are equivalent to a multi¬ 
dimensional radius in the configuration space, as a conversion to polar form will demonstrate. This 
transformation from Cartesian coordinates X to polar coordinates, i.e., the polar decomposition of 
the matrix X, was first executed in the context of a physical problem by Eckart. OTI and was first 
applied to polymer configuration problems by Solc. ll^ The transformation is the same as the 
well-known SVD decomposition in multivariate statistics. For our needs, the transformation of 
the volume element is most easily carried out with use of some matrix algebra. The Autonne- 
Eckart-Young theorem ll^ enables us to write X = R'^U where R G 50(3) is a rotation matrix 
in the special orthogonal (Eie) group {RR' = R'R = I3), ^ is a diagonal matrix such that is 
the matrix of eigenvalues of the non-singular matrix XX', and UU' = I3, U'U 7^ is a point in 
a Stiefel manifold. lIMll i.e., U lies in the coset space 0(/i)/0(/z — 3). The transformation of the 
3yU-dimensional volume element: dX —>■ J[dX/d{R^U)]dRd^dU, can be done in a few different 
ways. The Jacobian J[dX/d{R^U)] can be computed directly. |[3?1 but fewer computations are 
required if use is made of the rule from differential geometry that ^/gdx is the volume element 
associated to the metric = gijdxidxj . The Jacobian is computed in the Appendix. 

Eliminating the Zero Eigenvalue of the Laplacian 

Since the coordinates in eq. (7) are measured relative to the origin at the center of symmetry, the 
rows of X sum to zero. Define the orthogonal matrix T G 50(/i) that diagonalizes the Eaplacian 
via a similarity transformation. That is, make the substitution X = QT , such that tr(XKX') = 
tr(QTKT'Q') = tr(Qfi;Q'), where k is the matrix of non-vanishing eigenvalues of K, and Q is a 
normal coordinate representation of the configuration space. The zero eigenvalue of K is generated 
by the constant row of T , and the corresponding column in T' gives a zero when it multiplies 
X. This coordinate is the center of symmetry that is fixed at the origin, and XT' annihilates 
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the corresponding column of Q. One gets to the same place by deleting the constant column of 
T', such that T' —)■ T^; the zeros in Q and k combine and may be deleted, so that Q becomes 
a 3 X (/r — 1) matrix and k is a (/i — 1) x (/i — 1) positive definite diagonal matrix. Since 
XX' = QTqTqQ' = QQ' (note that TqTq = l^_i but TqTq ^ The Q coordinates can be 

reduced to polar form as described above for the cartesian matrix X, to give 

f3V = 7 /itr(SUfi;U') (9) 

where f/ is a 3 x (/r — 1) matrix representation of the coset space 0(/i — l)/0(/i — 4). The 
dependence on the rigid body rotation R cancels in eq. (9) because the trace is invariant to cyclic 
permutation of the arguments. 

The volume element on the configuration space begins as dX = 6{xo) where the 

delta function suppresses integration over the center of symmetry coordinate Xq. The Jacobian 
determinant of the transformation to normal coordinates is unity, and since the center of symmetry 
coordinate has been suppressed, we have 6{xo)dX ~ dQ. 

The transformation to polar coordinates (see Appendix) gives 

dX = dRdU[det{S)]^''-^^/^ n 1*^“ “ n 

a<b a 

to within an uninteresting constant. The integration over the space R G S'0(3) gives another 
constant, and since we are to integrate over dX/dS by the prescription of eq. (2), that is, the 
components of the diagonal tensor S are fixed, all that remains is the integral over U. Happily, the 
structure of k makes this part of the problem easier than one might expect. 

Evaluation of the Configuration Integral 

The probability that the James-Guth model elastomer is found in a state S = diag(S'a) is 

P{S)dS oc [det (^)](^-^)/2 TT 1 ^^ _ ^^ 1 ^^ f etr{-jfxSUKU')dU (10) 

a<b Juu'=i 

where etr(-) = exp[tr(-)]. Terms that would convert the proportionality to an equality will cancel, 
as does the free energy Aq, when the ratio P(S)(iS/P(So)(iSo = det(AA')P(S)/P(So) is com¬ 
puted. The integration over dU covers the coset space U G 0(/r — l)/0(/i — 4). The elements of 
U are bounded: 0 < |m“| < 1 because UU' = I 3 . (The extremely interesting geometrical fact that 
emerges from this coordinate transformation is this: In a configuration space of dN variables, three 
are deleted for center of mass motion, three for rigid body rotations, and three fix the macroscopic 
size of the system. The remaining dN — 9 variables comprise a non-Euclidean compact space. 
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independent of any model potential.) 

The evaluation of the integral is highly dependent on the eigenvalue speetrum, k, of the Lapla- 
eian. Several faets about the speetrum ean be dedueed from general theorems on matriees, but it 
is the small eigenvalues that are important. Regardless of the exaet values that the small eigenval¬ 
ues take, one knows that the dominant eontribution to the integral in eq. ( 10 ) will be attained in 
the region where UkU' has its minimum value. This means that the dominant eontribution to the 
integral will be found in the region of the U-space that is assoeiated with the smallest eigenvalues 
of K. In addition, the free energy must be an extensive thermodynamie funetion if the theory is to 
make any sense. But S = beeause of eq. (7). This requires that UkU' ~ 0(/i“^/^), so as 

to eaneel the maeroseopie size dependenee in S. These observations were made long ago;[l36l| the 
main purpose of the next seetion is to show that this estimate of the small eigenvalue is eorreet. 

Let the eigenvalues be ordered sueh that 0 < < K 2 < • • • < The minimum of 

ii{UKU') will be attained on the subspaee where U = {R, 0); R G 0(3) and 0 is a 3 X (/i 4) 

matrix of zeroes. In the next seetion we will see that the small eigenvalues are three fold degenerate, 
Ki = K 2 = ^ 3 , for a eubieal elastomer, whieh eollapses the integral to 

J etT{—jfiSUKU')dU ~ etr(— 7 /iKiS') = exp[— 7 /rKi(S'i + 52 + S 3 )] ( 11 ) 

This may be eonsidered to be the leading term in an expansion; it is diffieult to improve upon 
itlETl 

On evaluating the eontribution to the free energy from eq. (11) we get the term A = kT /iy^i (5i+ 
S2 + S3), where 7^1 is the only term that depends on eonstraints operating at the time of formation 
of the elastomer. These eonstraints are imprinted on the speetrum of eigenvalues of the Laplaeian 
K. The BFE average now eomes to the reseue: the average over the probability distribution of 
eonstraints is a simple average over eonneetivity. That is, (A) = kT /i( 7 Ki) {Si + ^2 + S3) requires 
a straightforward statistieal meehanieal average of the small eigenvalue. It is hard to imagine a 
teehnique for evaluating the small eigenvalues that does not invoke a BFE averaging proeedure. 
Reeognition that 7^1 is a BFE average, even with braekets omitted, is implieit in everything that 
follows. If the network is formed with a distribution of ehain lengths (MWD) the averaging of the 
small eigenvalue requires eonsiderable eare, and is left for another time. 

Estimating the Small Eigenvalues 

The . . smallest eigenvalue of the Laplaeian, is far from trivial; in fact, it is difficult to overem¬ 
phasize its importance.” {QjOt 

The smallest non-zero eigenvalue has been ealled the “algebraie eonneetivity” of a graph. [f38l 
Unfortunately, what is known about Ki in random graph theory literature is not very helpful; what 
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sets physical random graphs apart from mathematicians’ random graphs [|39l is that physical ran¬ 
dom graphs are embedding in with a more or less uniform density of vertices, and that has a 
profound influence on the spectrum at the small eigenvalue edge. For our purposes an embedded 
graph in a roughly isometric hard or soft solid has a number density p = p/V of vertices such 
that any slice through the graph with thickness I = and area A = 0(V^^^) contains 
vertices. For mathematical purposes this asserts a uniformity of density of vertices (to within the 
usual influence of the discontinuity at the surface) while allowing for short-range correlations that 
may be of interest (and will be of great interest for many problems). 

Using this notion of slices with thickness £, a cube with volume contains one vertex on 
average. The connections between a vertex and its neighbors depend, of course, on the detailed 
edge-length distribution in the graph. Suppose that a cubical elastic body is divided into a simple 
cubic lattice oi p = rA cells. A layer of these cells contains, on average, m? cross-links. The num¬ 
ber distribution of cross-links in every layer will be a Poisson distribution, with variance n, if the 
cross-links are distributed at random. Since we are interested in the large p behavior, the fluctua¬ 
tions in the number of vertices in the layers is sufficiently small to be neglected [variance/mean = 
0(?7,“^)]. Because these layers contain large numbers of vertices the statement effectively ignores 
short-range correlations. This is not an assertion that short-range correlations are unimportant. 

Our eigenvalue problem now inherits an enumeration scheme that looks like that used to cal¬ 
culate the vibrations in a simple cubic crystal. The elastomer does not have phonon modes in the 
sense of a crystalline solid. (However, the long wavelength modes are visible in a soft material as 
macroscopic oscillations or vibrations.) Instead, we have a connectivity problem that is mathemat¬ 
ically analogous to that of a crystalline solid. The long wavelength eigenvectors of the Laplacian 
matrix are the same as the long wavelength eigenvectors in a crystalline solid. It is well known 
that the long wavelength modes in a crystal determine the lowest energy excitations, and the long 
wavelength eigenvectors, with wavelength A S> max(£, (r^)^/^), will generate the small eigenval¬ 
ues of our graphs. James l[T3ll calculated some statistical properties of a regular cubic lattice in 
relation to his work on elasticity, so this idea is not new. What is new is a better way to handle the 
randomness that overrides the lattice. The fundamental reason for discretizing the space of the soft 
material is to borrow the matrix indexing scheme from crystal theory, where there is a one-to-one 
correspondence between coordinates and matrix indices; this greatly facilitates the construction of 
eigenvectors. There is no obvious indexing of a random graph, so the discretization is imposed to 
provide this correspondence. 

Index the cells in our subdivided cubical elastomer in the obvious way as j = {ji, j 2 , js}, where 
1 < ja < n. The probability p{r)dr ^ p(j“, that two average cells with multi-indices are 
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connected by a polymer chain is discretized as 


( 12 ) 

Summed over all eq. (12) yields the average number, /, of chains connected to the cross-link 
located at j". We are assuming a homogeneous elastomer, such that the number of chains con¬ 
nected to a cross-link selected at random is translationally invariant. In eq. (12) the parameter 
7 * = 3/2(r^)* is determined by the chain statistics at the time the network is formed. The fun¬ 
damental polymer physics that goes into the model elastomer asserts that (r^)* is the unperturbed 
mean-square end-to-end distance of a free chain at the temperature of cure. This is a fixed length 
parameter that will not vary with temperature. At the temperature of measurement, where we are 
computing the configuration integral in eq. (11), the corresponding parameter governing the po¬ 
tential of mean force is temperature dependent. At the time the network is formed the only role 
for the chain length distribution is to determine the connectivity. After the network is formed the 
chains assume their role of delivering stresses. 

The long wavelength (unnormalized) eigenvectors for a simple cubic lattice has components of 
the form 

exp[7rf(jimi -f ^ 2^2 -f 

where 1 < ja < 0 < rua < rUe, 1 < a < 3, and i = 1 /—1. The index set m = {mi, m 2 , m 3 } 

labels the eigenvalues The restriction me n limits the eigenvectors to long wavelengths; if 
the wavelength is less than several multiples of (r^)y^ the local structure of the graph will become 
important, and the corresponding eigenvectors will be different from simple Fourier functions (al¬ 
though one could write the exact eigenvectors as linear combinations of Fourier functions since the 
latter form a complete basis). The constant eigenvector, with mi = m 2 = m 3 = 0, generates the 
zero eigenvalue of K. The components of the eigenvectors for the simple cubic lattice are actually 
ria cos[(ja — l/2)ma7r/n]. The approximate complex version makes subsequent calculations sim¬ 
pler and does no damage. We have no need of toroidal boundary conditions, which would insert a 
factor of 2 in the trigonometric functions and artificially render the eigenvalues doubly degenerate. 

After going to the trouble of discretizing the space we will now undo this work by replacing 
sums by integrals. Of course, the motivation for discretizing is to make the physical picture clear 
and to help guide an understanding of the implications of the long wavelength regime where the 
method will be valid. Our eigenvector is now specialized to a particular plane wave for which 
mi = 1 and m 2 = m 3 = 0. We could just as well have chosen m 2 or m 3 to be the only non-zero 
index. That is, the isotropy of the network for a cubical elastomer renders the smallest non-zero 
eigenvalue triply degenerate. Since the selected eigenvector is constant on planes perpendicular to 
the x-axis, it picks up an entire layer of vertices in the slab of thickness and thereby averages the 
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connectivity over the slab. Our eigenvalue problem is now mapped into 


Kiti =Kti 

Ki exp{Titxi/ni) =f{exp{7nxi/ni) - j exp[- 7 ,(a;i - Xa)'] exp(7r^a;2/n£)da;2 (13) 

=/{l - (7,/7r) J exp[-74a;i - X 2 f] exp[- 7 rz(a;i - a;2)/n£]da;2} 

=/{l - exp[- 7 rV( 47 *n^^^)]} 

(There are cancelled factors of that arise from summing over the constant X 2 and x^ components 
of the eigenvectors on the left and over the cross-links in the slices perpendicular to the xi-axis on 
the right.) The first factor in brackets comes from diagonal elements of the matrix, and the second 
factor, the integral, arises from the average number / of off-diagonal elements, each of which has 
the value —1 in the Laplacian (the actual value will be —k, where k is the number of chains that 
connect the two cross-links in question to allow for multiple connectivity). 

Since n is as large as we like, it follows that the first several small eigenvalues will be given by 

= /{I - exp|-,r 2 mV( 47 ,n 2 (-")]} « 7 = 1117 ( 47 , 7 /’) = ,r=m=/{r=)./(6V'7’) (14) 

where = m\+ + m\] |m| < rUe . We now have the needed proof that the small eigenvalues 

are proportional to It is also clear that in the limit as k —)■ 0+, the spectral density g{K) 

of small eigenvalues of the Laplacian matrix tends to the same distribution that Debye calculated 
for a spherical continuum elastic solid. The only difference is that our long wavelength density is 
determined by microscopic parameters - the chemistry in / and chain statistics in 7 * - whereas 
Debye’s spectrum is determined by phenomenological Lame constants. Standard textbooks BOl 
explain the relation between the long wavelength spectrum for the Bom-von Karman model of a 
simple cubic lattice and the Debye low frequency spectrum. 

If there is a distribution of molecular weights (MWD) of the chains connecting cross-links the 
average over this distribution couples with the graph eigenvalue spectrum, requiring a very careful 
analysis of the averaging procedure. In first approximation this MWD average can be postponed 
to eq. (14), where it simply gives an average of 7 ^^ ~ 
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Summary 

On combining eqs. (1,3,11, and 14), and replacing the probability distribution by the equivalent 
Helmholtz free energy, one obtains 

A{S) =Ao(|^|) + [kTfifn\r^)^/{4{r^)/V^^/^)]{S, + S, + S^) 

- (/i - 5)kT\n{\S\^/^) - kT\n{\S, - ^211^2 - ^sH^i - ^sl) (15) 

Now that this result is at hand one can see that the last term in eq. (15), which causes the eigen¬ 
values of X to repel one another, is not of thermodynamic significance. It can contribute, say, a 
surface tension term only if at least one of the arguments — S';,! = 0[exp(/i“^/^)]. But this is 
ridiculously small; normal thermodynamic fluctuations in macroscopic lengths will be 0(/i“^/^), 
so that the term is at most 0[ln(/i)], and may be neglected. Furthermore, one can ignore the differ¬ 
ence between /i and /i — 5 in the preceeding logarithmic term. 

The next point to make is that the theory gives a stable un-stressed state. This state is a mini¬ 
mum of the free energy; that is 

aAisyas. = + g - ,kT/2s. = g - /^*=r/2-o.4o(|g|)/ai.,(|g|) 

d\S\ dSa Sa 

vanishes when evaluated at Sa = S^. Here G = A;T/r/7r^(r^)*/(4(r^)vf^^). At the minimum 

0 _ ^ikT/2 - {dAo{\S\)/dln{\S\)s=so _ fikT/2 + poV/2 

“ G G 

with the last version resulting from [S'! oc V'^. 

The term involving pq can be discarded with the following argument: Suppose we had a good 
theory for Aq, such that one could solve for the equilibrium volume {dAo/dV)T = —po = 0: the 
unstrained (uncompressed) state of the base polymer is a minimum when the hydrostatic pressure 
vanishes. The difference between the base polymer free energy and the approximate Aq that has 
been factored from the configuration integral in eq. (3) is owing to the constraints of the cross¬ 
links that are implicit in Aq, and this difference can be made small by decreasing the density of 
cross-links. The term with po may be set to zero, so that the only volume dependence that survives 
is that in G, so that = 0(1^^/^) as it must. What is most important about this observation is that 
the James-Guth theory does not have a physically realistic un-stressed state in the absence of the 
log(V^) term. 

The gyration tensor in the unstrained state of our cubical elastomer has equal diagonal compo- 
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nents 

rLo/2 

5 ° = V~^ / xldxidx2dx3 = Ll/12 

J—Lo/2 

giving Sa = A^(Lg/12) = /12. It also follows that |S'|/|S'°| = {V/Vy. Putting these 

pieees into eq.(15) gives 

AA = A/l„ + +\l + \l-3)- likT\n(,Ly (,6) 

The ratio eame from the eigenvalue speetrum and is fixed at the time of eure; this ratio 

is independent of temperature. Measurement of the stress-temperature eoeffieient gives d^) /dT, 
not dy)^,/dT. (This nomenelature differs from that used in other work so as to make elear that 
there are three states involved - the state of eure at T*, the unstressed state at temperature T, and 
the stressed state at temperature T'.) The number of ehains in the network is u = /i//2, and if one 
ignores the (usually very small) differenees between 14, V, andVo, the elastie free energy simplifies 
to 2 

= vkTiy/2L){^ )(-^i + ^2 + Ag — 3). 

Note that / differs from the maximum ehemieal funetionality so that it effeetively eorreets for 
dangling ehains and loops, and therefore v = /i//2 is the eount of “elastieally effeetive” ehains. 
The surprising presenee of tt in the modulus (front faetor) is a eonsequenee of the estimate of the 
small eigenvalue of the Laplaeian that was made with use of a plane wave eigenveetor. Of eourse, 
vr appears naturally in the Rouse speetrum of relaxation times; it should not be eonsidered unusual 
that it arises in the bulk elastomer eontext. Note that 7r^/24 0.41, whieh is is eomparable to 

values for this faetor that have appeared in the work by James and Guth and the later work of Flory 
invoking the eyele rank. 

Conclusion 

It is remarkable that the James-Guth many-body theory of elastieity gives essentially the same re¬ 
sult that Flory dedueed over a period of several deeades. The logarithmie term, ehain dimension 
ratio, and front faetor (7r^/24 rather than the eyele rank lfTTl l all emerge when the problem is han¬ 
dled in a natural way. As this work shows, what is natural requires some fairly intense mathematies 
(at least relative to the simpler Wall-Flory treatment), but the reward is a theory of elastieity that 
has a firm foundation. Very few pproximations have been made along the way; the most eritieal 
(and the hardest) plaee to make improvements is in the estimate of the small eigenvalue of the 
Laplaeian of the graph. The integral in eq. (11) might also be improved iHTll with deeper under¬ 
standing of the speetrum. Improvements in the potential alter the model, whieh goes beyond the 
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present concern. 

A final comment regarding potential improvements to the model potential needs to be made. If 
the many-body potential of mean force is to be augmented, for example with constraints, |[9l these 
must be inserted at the level of the coset space. That is, constraints depend on f/ — where 
U(g) is the locus of constraints, and not on X — X^. This is a critical observation because: (i) 
microscopic constraints cannot alter the macroscopic dimensions directly, although they have an 
indirect influence through the strength of microscopic forces, and (ii) the logarithmic term has to 
be maintained in the form of const, x /i In(l^) so as to yield a stable un-stressed state. 
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Appendix: Polar Coordinates for Matrices 

It seems that a heavy calculation is unavoidable to get the volume element for polar coordinates of 
matrices. The approach used here makes use of the well-known relation between the metric ds^ = 
Qijdxidxj, 1 < i,j < n for an n-dimensional real space and the volume element ^/gd:s., where 
dx = dxi. If coordinates are arranged in an m x n, m < n real matrix X the volume element 
Ilaj ^ < a < rn,l < j < n, and clearly this belongs to the metric ds‘^ = ii{dXdX') . 
(Note that the rank of X is assumed to be m throughout this discussion.) The differential of the 
polar decomposition. If^ [34l . X = R!^U, is 

dX = R'{6r^U + diU + ^dU) (17) 

where 6 r = RdR' is a skew-symmetric matrix. Here R G 0{m) spans the group of orthogonal 
matrices, are the eigenvalues of XX', and UU' = Im, but U'U 7 ^ if m 7 ^ n; in the latter 
case (/ is a Stiefel manifold, l.e., the coset space U E 0{n)/0{n — m), and if m = n, (7 G 0(n). 
Inserting eq. (17) into the metric and simplifying gives 

= tr[d.^^ + {6r^ + ^6u){^6r' + 5u'^) + ^‘^dU{l — U'U)dU'] (18) 

Similar to 5r, 5u = UdU' is skew-symmetric. If m = n the last term vanishes. 

The middle term in eq. (18) is simplified by extracting symmetric and skew-symmetric parts 
with 5rS,+i5u = (.^d^-f -f {^56 — 59'^), where dcf) = (dw — dr)/2 and 59 = (dw -|-dr)/2 are 
both skew-symmetric. But now the first term in parentheses is symmetric and the second is skew- 
symmetric, which enables one to simplify the expression tr [(s -f a) {s' -f a')] = tr [(s -|- a) (s — a)] = 
tr(s^ — a^) where s and a are the symmetric and skew-symmetric, respectively. In component 
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form, this term becomes 


tr((ir,^ + ^6u){^6r' + 6u'^) — 2 — ibf‘^(t>^ab + ('Ca + 

a<h 

with associated volume element 

I'Ca {d^pabdOab ~ I'Ca \ dTahdUah-i 

a<b a<b 

where numerical factors of no consequence to thermodynamics have been dropped. 

The last term in eq. (18) requires some work, but is instructive. Write a partitioned U as 
U = {Ui,U 2 ) = Ui{l,y), where y = f/f^f/ 2 - The matrix Uiis m x m and is non-singular 
except on subspaces of lower dimension. Since UU' = 1 it follows that UU' = f/i(l + yy')U[ —)■ 
1 + yy' = {U[Ui)~^. Some algebra establishes 

1 - j (1 + y'yr\y\ -i) = (y, -i)c(i + y'y)-\y\ -i) 

with use of y'{l + yy') = (1 + y'y)y' and its variations. The dimension of the ^/-manifold, a 
Grassmannian, is m(n—m). Since df/ = (df/i, df/iy+f/id?/) it follows that —l)c = —Uidy, 
so that 

tr| [i^^dUil - U'U)dU'] = tr| [U'^i^Uidy{l + y'y)~^dy'] (19) 

The final detail that one needs to calculate the volume elements is this: Let dZ = (dzaj), 1 < a < 
m, 1 < j < n, be an m X n matrix with row form 

dz = {dzii, dzi2, ■■■ , dzin, dZ2l, ■■■ , dZ2n, ■■■ , dZmn)] 

it can be seen that tr{AdZBdZ') = dz{A' G) B)dz', where A and B are conformable. The volume 
element associated with this metric is dV = dzaj. The volume element associated 

to the component defined in eq. (19) is 

|t/;{2c/,|(»-™)/2|l + y'y\-'^l^\'[\dvj =|«2|(»—)/2|l + + !/Vr’"''TI‘'!'« 

=l?¥”-’”)/"|i + ro'r”/Tl‘*!'“^ 
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Putting the pieces together one finds 


JJdxaj =const.{|^^|(” ib\Wdia}{drdu}{\l + yy'\ JJ JJ d|/afc 

a,j a<b a a=l k=l 

=const.{|e"|("-”^)/2 JJ 1^2 _ ^2\Yld^a}dRdU (20) 

a<b a 

where dR and dU are short-hand notations for the volume elements on SO{m) and SO{n) / SO{n— 
m), respectively. A check of physical dimensions (L = length) will show that dim((iX) = 

a dia\ ~ since the orthogonal group and the Steifel manifold 

are dimensionless. 
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